{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation problems\n",
    "\n",
    "In this tutorial, we explain one of the greatest features of GeoStats.jl: the ability to setup geostatistical problems indenpendently of the solution strategy.\n",
    "\n",
    "If you are an experienced user of geostatistics or if you do research in the field, you know how hard it is to compare algorithms fairly. Often a new algorithm is proposed in the literature, and yet the task of comparing it with the state of the art is quite demanding. Even when a comparison is made by the author after a great amount of effort, it is inevitably biased.\n",
    "\n",
    "Part of this issue is attributed to the fact that a general definition of the problem is missing. What is it that we call an \"estimation problem\" in geostatistics? The answer to this question is given below in the form of code.\n",
    "\n",
    "Before we proceed, please install the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStats is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage Plots is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage PyPlot is already installed\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "for pkg in [\"GeoStats\", \"Plots\", \"PyPlot\"]\n",
    "    Pkg.add(pkg)\n",
    "end\n",
    "\n",
    "# make sure this tutorial is reproducible\n",
    "srand(2000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem definition\n",
    "\n",
    "An estimation problem in geostatistics is a triplet:\n",
    "\n",
    "1. Spatial data (i.e. data with coordinates)\n",
    "2. Spatial domain (e.g. regular grid, unstructured grid)\n",
    "3. Target variables (or variables to be estimated)\n",
    "\n",
    "It doesn't involve variograms, training images, or any tuning parameter. These concepts belong to solvers.\n",
    "\n",
    "Let's make it a concrete example, and create some data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "3×4 GeoDataFrame (x and y)\n",
       "<table class=\"data-frame\"><thead><tr><th></th><th>x</th><th>y</th><th>station</th><th>precipitation</th></tr></thead><tbody><tr><th>1</th><td>25.0</td><td>25.0</td><td>palo alto</td><td>1.0</td></tr><tr><th>2</th><td>50.0</td><td>75.0</td><td>redwood city</td><td>0.0</td></tr><tr><th>3</th><td>75.0</td><td>50.0</td><td>mountain view</td><td>1.0</td></tr></tbody></table>"
      ],
      "text/plain": [
       "3×4 GeoDataFrame (x and y)\n",
       "\n",
       "│ Row │ x    │ y    │ station       │ precipitation │\n",
       "├─────┼──────┼──────┼───────────────┼───────────────┤\n",
       "│ 1   │ 25.0 │ 25.0 │ palo alto     │ 1.0           │\n",
       "│ 2   │ 50.0 │ 75.0 │ redwood city  │ 0.0           │\n",
       "│ 3   │ 75.0 │ 50.0 │ mountain view │ 1.0           │"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStats\n",
    "\n",
    "geodata = readtable(\"data/precipitation.csv\", coordnames=[:x,:y])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `GeoDataFrame` type is defined in GeoStats.jl as a lightweight wrapper over Julia's DataFrame types. Besides the data, it stores the name of the columns representing spatial coordinates. It is also possible to load the data directly into a `GeoDataFrame` by using the `readtable` function, and passing the option `coordnames=[:x,:y]`.\n",
    "\n",
    "Next, we define the domain in which the variables will be estimated. One of the many choices possible is the regular grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RegularGrid{Float64,2}\n",
       "  dimensions: (100, 100)\n",
       "  origin:     (0.0, 0.0)\n",
       "  spacing:    (1.0, 1.0)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "domain = RegularGrid{Float64}(100,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that by default, a regular grid has zero origin and unit spacing. Alternatively, these parameters can be specified explicitly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RegularGrid{Float64,2}\n",
       "  dimensions: (100, 100)\n",
       "  origin:     (0.0, 0.0)\n",
       "  spacing:    (1.0, 1.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "domain = RegularGrid([100,100], [0.,0.], [1.,1.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like many other domain types in GeoStats.jl, regular grids are lightweight objects. They do not require any memory space other than the space used to save the input parameters (i.e. dimensions, origin and spacing). Looping over a regular grid or finding the coordinates of a given location is done analytically for maximum performance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we define the estimation problem for the precipitation variable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D EstimationProblem\n",
       "  data:      3×4 GeoDataFrame (x and y)\n",
       "  domain:    100×100 RegularGrid{Float64,2}\n",
       "  variables: precipitation (Float64)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "problem = EstimationProblem(geodata, domain, :precipitation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the problem\n",
    "\n",
    "Now that the problem is unambiguously defined, we can solve it with various estimation solvers. In this tutorial, we will use the polyalgorithm Kriging solver distributed with GeoStats.jl. In Kriging, each variable of the problem is (optionally) parametrized by a mean and a variogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GeoStats.Kriging\n",
       "  └─precipitation\n",
       "    └─variogram = GeoStats.GaussianVariogram{Float64,Float64,Distances.Euclidean}\n",
       "  range: Float64 35.0\n",
       "  sill: Float64 1.0\n",
       "  nugget: Float64 0.0\n",
       "  distance: Distances.Euclidean\n",
       "\n",
       "    └─mean = nothing\n",
       "    └─degree = nothing\n",
       "    └─drifts = nothing\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solver = Kriging(\n",
    "    :precipitation => @NT(variogram=GaussianVariogram(range=35.))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The line above translates to \"solve the precipitation variable using a Gaussian variogram\". The prefix `@NT` extends for `NamedTuple` and won't be necessary in Julia v0.7 and beyond. When only the variogram is specified, Ordinary Kriging is triggered.\n",
    "\n",
    "The user can specify the mean (e.g. `mean=.5`) for Simple Kriging, the polynomial degree (e.g. `degree=1`) for Universal Kriging, and the drift functions (e.g. `drifts=[x -> 1 + x[1], x -> 2x[2]]`) for External Drift Kriging. For more solver options, please consult the GeoStats.jl documentation.\n",
    "\n",
    "The solution to the problem is easily obtained with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D EstimationSolution\n",
       "  domain: 100×100 RegularGrid{Float64,2}\n",
       "  variables: precipitation"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = solve(problem, solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is stored in an efficient format with all the necessary information to reconstruct the estimates spatially. The user can digest the solution into a Julia dictionary by calling `digest(solution)`. The resulting dictionary maps variable names (e.g. :precipitation) to their mean and variance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{Symbol,Dict{Symbol,Array{Float64,2}}} with 1 entry:\n",
       "  :precipitation => Dict(:mean=>[0.687121 0.689082 … 0.671328 0.671373; 0.68908…"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = digest(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, very often we just want to visualize the results. In GeoStats.jl, solutions to estimation problems can be plotted directly in a standardized format for comparison between different solvers and parameter settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots; pyplot()\n",
    "\n",
    "plot(solution, size=(1000,500))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "The ability to work at the level of the problem definition is quite desirable. Users can switch between different solvers without having to learn new syntax. This approach also guarantees that the problem being solved is the same and that the comparison is fair.\n",
    "\n",
    "Documentation on how to write estimation solvers can be found in the [GeoStatsBase.jl](https://github.com/juliohm/GeoStatsBase.jl) package. Additional helper tools for developers can be found in [GeoStatsDevTools.jl](https://github.com/juliohm/GeoStatsDevTools.jl). Please open an issue in any of these modules if you have questions about solver development."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
